#include "gc.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "Path.h"

typedef std::list<Point> PointList;
typedef std::vector<PointList> PointListSet;

void evaluate (std::string truthPathName, int truthPathNo, 
	       std::string truthImageFileTemplate, 
	       std::string resPathName, int resPathNo, 
	       std::string resImageFileTemplate)
{
  ListOfPaths truthPathList, resPathList;
  truthPathList.readXMLPaths(truthPathName.c_str());
  resPathList.readXMLPaths(resPathName.c_str());
  Path truthPath = truthPathList.get(truthPathNo);
  Path resPath = resPathList.get(resPathNo);
  int truthStartSliceNo = truthPath.lastSlice();
  int truthEndSliceNo = truthPath.firstSlice();
  int truthSliceNum = truthEndSliceNo - truthStartSliceNo + 1;
  int resStartSliceNo = resPath.lastSlice();
  int resEndSliceNo = resPath.firstSlice();
  int resSliceNum = resEndSliceNo - resStartSliceNo + 1;
  PointListSet truthPointListSet(truthSliceNum);
  PointListSet resPointListSet(resSliceNum);
  for (int truthIndex = truthPath.size() - 1; truthIndex >= 0; truthIndex--) {
    int sliceNo = truthPath.getSliceAtIndex(truthIndex);
    int label = truthPath.getLabelAtIndex(truthIndex);
    char fn[200];
    sprintf(fn, truthImageFileTemplate.c_str(), sliceNo);
    LabelImage::Pointer image = readLabelImage(fn);
    itk::ImageRegionConstIteratorWithIndex<LabelImage> 
      itr(image, image->GetLargestPossibleRegion());
    while (!itr.IsAtEnd()) {
      if (itr.Get() == label) {
	LabelImage::IndexType index = itr.GetIndex();
	Point point(index);
	truthPointListSet[sliceNo - truthStartSliceNo].push_back(point);
      }
      ++itr;
    }
  }
  for (int resIndex = resPath.size() - 1; resIndex >= 0; resIndex--) {
    int sliceNo = resPath.getSliceAtIndex(resIndex);
    int label = resPath.getLabelAtIndex(resIndex);
    char fn[200];
    sprintf(fn, resImageFileTemplate.c_str(), sliceNo);
    LabelImage::Pointer image = readLabelImage(fn);
    itk::ImageRegionConstIteratorWithIndex<LabelImage> 
      itr(image, image->GetLargestPossibleRegion());
    while (!itr.IsAtEnd()) {
      if (itr.Get() == label) {
	LabelImage::IndexType index = itr.GetIndex();
	Point point(index);
	resPointListSet[sliceNo - resStartSliceNo].push_back(point);
      }
      ++itr;
    }
  }
  // // For debug
  // std::cout << "Truth:" << std::endl;
  // for (int sliceNo = truthStartSliceNo; sliceNo <= truthEndSliceNo; sliceNo++) {
  //   std::cout << "Slice " << sliceNo << ": " 
  // 	      << truthPointListSet[sliceNo - truthStartSliceNo].size()
  // 	      << std::endl;
  // }
  // std::cout << "Res:" << std::endl;
  // for (int sliceNo = resStartSliceNo; sliceNo <= resEndSliceNo; sliceNo++) {
  //   std::cout << "Slice " << sliceNo << ": " 
  // 	      << resPointListSet[sliceNo - resStartSliceNo].size()
  // 	      << std::endl;
  // }
  // // ~ For debug
  int truthSliceNo = truthStartSliceNo;
  int resSliceNo = resStartSliceNo;
  int tp = 0;
  int fp = 0;
  int fn = 0;
  while (truthSliceNo <= truthEndSliceNo && resSliceNo <= resEndSliceNo) {
    while (!truthPointListSet[truthSliceNo - truthStartSliceNo].empty()) {
      Point truthPoint = truthPointListSet[truthSliceNo - truthStartSliceNo].front();
      truthPointListSet[truthSliceNo - truthStartSliceNo].pop_front();
      bool isFound = false;
      for (PointList::iterator itr = 
	     resPointListSet[resSliceNo - resStartSliceNo].begin(); 
	   itr != resPointListSet[resSliceNo - resStartSliceNo].end(); itr++) {
	if (itr->isEqual(truthPoint)) {
	  tp++;
	  resPointListSet[resSliceNo - resStartSliceNo].erase(itr);
	  isFound = true;
	  break;
	}
      }
      if (!isFound) {
	fn++;
      }      
    }
    fp += resPointListSet[resSliceNo - resStartSliceNo].size();
    truthSliceNo++;
    resSliceNo++;
  }
  for (; truthSliceNo <= truthEndSliceNo; truthSliceNo++) {
    fn += truthPointListSet[truthSliceNo - truthStartSliceNo].size();
  }
  for (; resSliceNo <= resEndSliceNo; resSliceNo++) {
    fp += resPointListSet[resSliceNo - resStartSliceNo].size();
  }
  // std::cout << "TP = " << tp << std::endl;
  // std::cout << "FN = " << fn << std::endl;
  // std::cout << "FP = " << fp << std::endl;
  std::cout << "Precision = " << double(tp) / double(tp + fp) << std::endl;
  std::cout << "Recall = " << double(tp) / double(tp + fn) << std::endl;
} 



int main (int argc, char* argv[])
{
  int argNo = 1;
  std::string truthPathName = argv[argNo++];
  int truthPathNo = atoi(argv[argNo++]);
  std::string truthImageFileTemplate = argv[argNo++];
  std::string resPathName = argv[argNo++];
  int resPathNo = atoi(argv[argNo++]);
  std::string resImageFileTemplate = argv[argNo++];
  evaluate(truthPathName, truthPathNo, truthImageFileTemplate, 
	   resPathName, resPathNo, resImageFileTemplate);
  return 0;
}
